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Abstract x s liquid height above the jet 

The paper presents a numerical solution of an y distance measured from the nozzle wall 

axisymmetric turbulent jet discharging axially from 
below into a cylindrical tank, and directed towards Greek symbols: 
the liquid vapor interface. The liquid-vapor inter- 
face is assumed to be flat and shear free. The e dissipation rate of turbulent kinetic 

k-e turbulence model is used to calculate the eddy energy 

viscosity. The turbulence i ntensi ty di stri bution 

and the length scale associated with the k-e model k Von Karman constant (= 0.42) 

are calculated as functions of jet flow rates and 

system parameters. Numerical results are compared A integral length scale of turbulent 

with appropriate experimental data. The problems eddies 

associated with the free surface boundary condi- 
tions for the turbulent quantities are discussed. u liquid viscosity 

Nomenclature eddy viscosity 

Ci, C 2 , constants in k-e turbulence model Ueff effective viscosity 

D tank diameter p liquid density 

d jet or nozzle diameter Prandtl number in k equation 

G turbulence generation term a e Prandtl number in e equation 

g acceleration due to gravity Subscripts: 


k turbulent kinetic energy 


evaluated based on jet nozzle dimension 


L 

1 
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0 

R 
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U 

u 
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length scale, k^ 2 /e 
Prandtl mixing length 
pressure 

volume flow rate 
tank radius 

radial coordinate measured from the 
centerline 

jet Reynolds number, pUjd/p 
mean axial velocity 

rms value of axial fluctuation velocity 

mean radial velocity 

rms value of radial fluctuation 
velocity 

axial coordinate measured from the 
tank bottom 


s evaluated at liquid-vapor interface 

location 


Background 


The present work is motivated by the fluid mix- 
ing and the liquid-vapor interface condensation 
problems which may be encountered during the on- 
orbit storage and transfer of subcritical cryogens 
and the pressure control of the storage tank. The 
liquid-vapor interface location and configuration 
are not well defined in micro-gravity environments. 
Therefore, direct venting of vapor to control the 
tank pressure either during long-term storage or 
during the filling process of an empty tank is not 
a viable option. The problems associated with 
direct venting in low-gravity conditions are ade- 
quately described by Labus et al. 1 Realization of 
these problems has led to the development of the 
concept of thermodynamic vent system 2-5 (TVS) for 
the pressure control of cryogenic storage tanks and 
the "no-vent" fill technique 5 ’** for filling tanks 
with cryogens. The rate of vapor condensation on 
to a liquid surface and the fluid mixing process 
are important problems associated with these 
systems. 


K »' . 
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Fluid Mi xinq 

Experimental investigations 7-9 of fluid mixing 
due to an axial dyed jet introduced into an isother- 
mal fluid (ethanol), were performed in both normal 
and reduced gravity environments. Based on the 
findings of these tests, axial liquid jets are con- 
sidered as the preferred technique for providing 
fluid mixing in low-gravity conditions. The 
desired objective of the use of axial liquid jet is 
to control the tank pressure by destroying thermal 
stratification of the stored liquid and by enhancing 
the rate of condensation of the vapor on the liquid- 
vapor interface. 

Vapor Condensation on the Liquid Surface 


Liquid is withdrawn from the outer part of the bot- 
tom of the tank. The volume flow rate of the 
liquid withdrawn is the same as that of the injected 
liquid. Therefore, the liquid level above the jet 
remains constant. The outflow area is made much 
larger than the inflow jet area, so that the outflow 
velocity is negligibly small compared with the jet 
velocity. However, in numerical computation the 
outflow velocity term is retained, even though it 
is two orders of magnitude lower than the jet inlet 
velocity. The injected liquid and the tank liquid 
temperatures are assumed to be the same and uniform. 
Therefore, the effects of heat transfer and thermal 
stratification on the flow field and the evolution 
of turbulence are not considered in the present 
analysis . 


Experimental findings have established that the 
rate of the vapor condensation onto the liquid can 
be greatly enhanced by the introduction of turbu- 
lence beneath the liquid surface. Sonin et al. 10 
and Helmick et al. 11 measured the steady state rate 
of condensation of nearly saturated steam onto a 
turbulent liquid surface. The turbulence on the 
liquid side was produced by a single axial jet 
directed towards the liquid surface from below. 

The geometrical configuration shown in Fig. 1 is an 
adaptation of the test cell used by Sonin et al. 10 
For a fixed water level, condensation rates were 
experimentally determined as a function of liquid 
side turbulence intensity. The turbulent intensity 
was controlled by varying the jet flow rate. The 
turbulence intensity v (rms value of fluctuating 
velocity), its spatial distribution and the turbu- 
lence macro scale, were measured in the absence of 
condensation. These measurements were performed at 
axial locations greater than 3D from the jet inlet. 
It was argued and subsequently supported by experi- 
mental measurement, that in the upper region (x/D > 
3), the turbulent velocities are much larger than 
any mean circulatory velocities. It was also estab- 
lished that at some distance below the interface 
(away from the damped layer) the turbulence is 
approximately isotropic in the horizontal plane and 
is essentially unaffected by the proximity of the 
interface, and the length scale of turbulence is 
locked to the system diameter, D. The most signifi- 
cant finding of these investigations 10 * 11 is that 
for a shear-free and wave-free liquid interface, 
the condensation rate can be characterized in terms 
of the liquid-side turbulence characteristics v 
and A. Empirical correlations to describe the con- 
densation rate, the turbulence intensity near the 
surface, and the conditions under which they apply 
are fully discussed in Refs. 10 and 11. The objec- 
tive of the present work is to compute turbulence 
characteristics of the physical system of Refs. 10 
and 11 by solving the appropriate governing equa- 
tions numerically. 

Problem Formulation 


The Physical Problem 

The geometric configuration and the coordinates 
system used to analyze the problem are shown in 
Fig. 1. A cylindrical tank of diameter D is par- 
tially filled with liquid. A circular liquid jet 
of diameter d directed upward along the tank cen- 
tral axis is introduced at the bottom of the tank. 


The geometric parameters of the problem are 
the ratios of the tank diameter to jet diameter, 

D/d, and the liquid height above the jet to the 
tank diameter x s /D. In the present analysis only 
the parameter x s /D is varied. 

Governing Equations 

The jet considered in the present analysis is 
an axi symmetri c turbulent jet. The flow is steady- 
state and incompressible. The liquid-vapor inter- 
face is assumed to be flat and shear free. The 
k-e model of Jones and Launder 12 is chosen for cal- 
culating the eddy viscosity. The time-averaged con- 
tinuity, momentum and k-e equations are: 

Conservation of mass 

fj <pU> + ^ <r p v> - 0 (1) 

Conservation of x-momentum 

fe K> * 4 M • fe Gv, I) * Hr (%rr » 
* {- 8 * « * fe Ovr S - ! <*) * H? ( r >‘eff 1)1 


Conservation of r-momentum 


fe H * Hr ( kp » 2 ) - 1; (v™ I) 

•HrOvrfiM-lf-fcGwB 


( ril eff 3r " I prk ) ‘ Zv eff 'l | <3) 


T r9r 

k-e turbulence model 


fe <»«> . H? < pp «> - fe 


p + 


M 3k 


r3r 




HU 3k 
a. 


+ {p eff G - pc) 


(4) 
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fc <pUe> * tJf (rpUO 


3_ 

ax 



M 3e 
aj 9x 



where the turbulence production term, G, is defined 
by 



The terms in the braces, on the right-hand side of 
each equation are the source terms. The turbulent 
viscosity is given by the expression: 

•*t = Si fi r ; M eff - v + (7) 

Numerical values of the empirical constants associ- 
ated with the turbulence model are 12 * 13 

C ] -1.44, 

°k = K0 ' 

( 8 ) 


C 2 -1.92 , 


C = 0.09 , 


k * 0.42 and o 


’ (°j - c .) 


. 1/2 


Boundary Conditions 

The flow domain between 0 < r < 0/2, and 
0 < x < x s is considered for the computations. 

The flow problem considered in the present analysis 
has a line on symmetry along the tank central 
axis. The symmetric conditions at the center line 
are: 

3U 3k 3e n , v/ A /cn 

3? “ 9? = IF = 0 * and V “ 0 (9> 

At the wall no-slip boundary conditions apply 

U - v - o . <10) 


Also, the wall functions 14 are used to calculate 
shear- stress and are included as a source term for 
velocity component parallel to the wall in the 
governing equations. 15 

The velocity and the turbulent kinetic energy 
are assumed uniform at the jet inlet. 

U - U. , V - 0 and k - 0.005 U? (11) 

J J 

The distribution of the dissipation rate is given 
by 6 : 


e 


c 3/4 k 


3/2 

~ 


( 12 ) 


Where the mixing length scale, 1, is given by the 
expressions: 
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1 = icy for y < where y = | - r 

1 ■ c ,(!) f ° r » i (£>(/-) ■ <13> 

At the bottom of the tank (i.e., x = 0) , where 
the liquid is withdrawn, the velocity is assumed to 
be uniform and equal to Q/A, where Q is the jet 
inflow rate and A is the tank outflow area. 

The boundary conditions for the turbulence 
quantities, k and e at the liquid vapor inter- 
face, pose special problems in the numerical compu- 
tation. To our knowledge, the appropriate boundary 
conditions for k and e at the interface have 
not been established either theoretically or empiri- 
cally. In the present analysis it is assumed that 
at the free surface 


3k 

3x 



(14) 


Also, the assumption of the flat and shear free 
interface implies: 

U - 0 . = 0 at x=x s (15) 


Numerical Solution 


Solution Methodology 

The governing equations (Eqs. (1) to (5)) were 
solved by a finite difference numerical scheme. 15 ’ 15 
The finite-difference equations were derived by 
integrating the differential equations over an ele- 
mentary control volume surrounding a grid node 
appropriate for each dependent variable. A stag- 
gered grid system was used. The scalar properties 
P, k and e were stored midway between the U 
and V velocity grid nodes. The bounded skew 
hybrid differencing 16 was incorporated for the con- 
vective terms and the integrated source terms were 
linearized. The pressure distribution was obtained 
by means of pressure implicit split operation 
(PISO) 17 predictor-corrector technique. 

The calculations were performed with a nonuni- 
form grid distribution with the concentration of 
grid nodes near the center line, the wall and the 
interface regions where the gradients of flow prop- 
erties are expected to be large. The grid independ- 
ent solutions were determined by trial runs. 

A 36 by 31 grid in the x and r direction, 
respectively, was used for x s /D = 3.67. Comparison 
of numerical results of this grid size with that 
obtained by using higher grid numbers shows that 
the solutions are essentially grid independent. 

For x s /D greater that 3.67, higher grid numbers 
were used in the x-direction (56 by 31, for 
x s /D = 7.2). 

The finite difference equations were solved 
iteratively by the ADI method with under relaxa- 
tion 15 until the solution converged. The maximum 
absolute residual sums for each dependent variable 
less than 0.005 was used as a criterion of conver- 
gence. Based on this criterion, the convergent 


3 


solutions were obtained within 1500 iterations, 
with 0.8 sec/iteration on CRAY-XMP CPU time for 
initial guessed solutions set of zero. 

Parameters 

The numerical solutions presented in this 
paper were performed using the input conditions 
representative to those of experiments of Refs. 10 
and 11. The thermophysical properties of liquid 
water at 20 °C were used in all computations. The 
tank diameter, D, and the jet diameter, d, were set 
equal to 15.3 and 0.64 cm, respectively (same as 
those of Ref. 10). The jet Reynolds numbers were 
varied between 30 000 to 110 000, depending upon 
the liquid height to the tank diameter ratios, 
x s /D. Most of the calculations were performed for 
x s /D of 3.67 and 4.65. For high jet Reynolds 
number, numerical solutions were also obtained for 
x $ /D = 7.2, to investigate what effect, if any, the 
large liquid depth has on the solution. 

Results and Discussion 

The numerical solutions of the differential 
Eqs. (1) to (5), with the empirical coefficients of 
the turbulence model given by Eq. (8), subject to 
the boundary conditions (9) to (15), yielded distri- 
butions of mean flow velocities U and V, the tur- 
bulent kinetic energy, k, and the dissipation of 
the turbulence kinetic energy, e. As noted in the 
preceding section, the geometrical configuration and 
the input conditions used in the present analysis 
are similar to those of the experiments. 10 *^ 
Therefore, the numerical results that offer direct 
comparison with the experiments are presented in 
this paper. A detailed analysis of the numerical 
results encompassing a wide range of parameters is 
still being performed and will be reported in the 
future. 


Figure 3 shows the variation of the turbulence 
velocity along the central axis as a function of 
distance from the interface. The measured values 
of the vertical and horizontal fluctuating veloci- 
ties as reported in Ref. 11 are included for 
comparison. Both the measured and calculated 
results correspond to jet Reynolds number of 51 000 
and x s /D * 3.67. The experimental results suggest 
that a depth of about 0.13D below the interface 
the turbulence is approximately isotropic and is 
essentially unaffected by the interface. The numer- 
ical results are in excellent agreement with the 
experiment in this region. The vertical turbulent 
velocity closer to the surface damps out. Its 
kinetic energy is imparted to the horizontal compo- 
nents with a resultant increase in the horizontal 
turbulent velocities. Therefore, the turbulence 
closer to the surface is nonisotropic and the use 
of k-e model may not be appropriate. The compari- 
son between the numerical results and experimental 
data worsens as the free surface is approached. 

Sonin, et al. 10 obtained an approximate analyt- 
ical solution for the turbulence distribution by 
applying the high Reynolds number form of the k-e 
model of turbulence to the region x/D > 3, where 
the mean circulatory velocity is small compared 
with fluctuating velocities. They assumed that the 
turbulent kinetic energy, k, can be approximated as 
a function of x alone and the length scale 
L = k 3/2/ e in the model is constant and locked to 
the tube diameter. The analytical solution of the 
k-e model for purely diffusive turbulence shows 
that the isotropic turbulent velocity and the 
length scale are given by the expressions: 

* ■ '(LHif) <l6 ’ 

and 


The numerical results for the rms value of the 
turbulence velocity, v, (calculated from the expres- 
sion k = (3/2)v 2 ) at a fixed axial location, 
x/D = 3.33, and two radial locations r/R = 0 and 
r/R = 0.62, are presented in Fig. 2, where they are 
plotted as a function of the jet flow rate. The 
abscissa is chosen as Q/Dd, to facilitate compari- 
son with the experiment. The experimental rms 
values 10 of either the vertical or horizontal fluc- 
tuating velocity are shown on the figure for 
comparison. The range of jet Reynolds numbers for 
the stated range of flow rates varies between 
30 000 to 80 000. The liquid height to the tank 
diameter ratio, x s /D, for the experiment was 3.67. 
The numerical results presented in Fig. 2 correspond 
to this value of x s /D. The horizontal sampling 
region for the experimental data extended from the 
axis to the 2/3 radius of the test cell. The exper- 
imental data show that for x/D > 3, and at a 
istance of the damped layer thickness (*0.3D) below 
the interface, the turbulence is essentially iso- 
tropic in a horizontal plane. The numerical results 
obtained by using the k-e model are shown in 
Fig. 2. The solid lines drawn in Fig. 2 show the 
numerical results of v at the axis and at a radial 
distance of 0.62R. Almost all of the measured data 
fall between these two lines and show excellent com- 
parison with the numerical results. Both the exper- 
imental as well as the numerical results show almost 
a linear relation between v and Q/Dd. 


k 3/ 2 

L = = |3D where b and |3 are constants 

e 


(17) 

By comparing this solution with the experimental 
results Sonin et al. 10 suggested the values of 
13-1.1 and b = 21.8 for high Reynolds number 
flows . 

The experimental data 10 used to fix the value 
of (3 are reproduced in Fig. 4 along with the 
numerical solution of the corresponding problem. 

The numerical results show that the slope of 
ln(v) versus x/D is not constant, which implies 
that 6 is not a constant in the range of 
3.1 < x/D < 4.2 as suggested in Ref. 10. Analysis 
of numerical results shows that for x/D, up to 
approximately 3.6, the generation term in k-e 
equations is of the same order of magnitude as the 
diffusion term. Therefore, the solution obtained 
for purely diffusive turbulence may not apply in 
this region. The length scale, L is plotted in 
Fig. 5 as a function of the radial distance from 
the axis. Figure 5 shows that the length scale at 
a given axial location is essentially uniform 
except close to the tank wall. However, the length 
scale varies with the axial distance which is not 
in agreement with Sonin, et al. 10 conclusion that 


I 
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for high Reynolds number the length scale is con- 
stant and locked to the system diameter. Near the 
surface, the numerical solution for the jet Reynolds 
number of 50 000 shows that L/D varies between 

0.42 and 0.6. 

The analytical solution obtained in Ref. 10 
should indeed be valid if x is sufficiently large 
(theoretically x -► . In order to investigate 

the effects of the liquid depth and the boundary 
conditions (14) for k and e imposed at the 
interface, numerical solutions for a fixed jet 
Reynolds number of 110 000 were obtained for liquid 
depths of x s /D = 4.65 and 7.2. The variation of 
the turbulent kinetic energy, k, and the length 
scale, L, as a function of x/D is shown in 
Figs. 6 and 7. As seen from these figures the solu- 
tion remains practically unaffected for x/D up to 
3.6, and the length scale L/D attains almost a 
constant value 0.9 for large x/D. The region 
3.6 < x/D < 4.6 probably represents transition to. 
purely diffusive turbulence where the effect of mean 
circulatory is negligible. However, the numerical 
results show that turbulent kinetic energy in steady 
diffusive turbulence region increases or approxi- 
mately remain constant with the axial distance, x, 
which is unrealistic. In our opinion, it is caused 
by the imposition of the boundary condi tons of 3/8x 
of all turbulence quantities zero at the interface. 
Also the constants incorporated in the model may not 
be appropriate. In the damped layer, close to the 
surface where the vertical component of turbulent 
velocity damps out, the use of the Reynol ds-Stress 
equations may be useful. However, before a particu- 
lar model is chosen, we feel that it is more impor- 
tant to resolve first, the issues associated with 
the free surface boundary condtions for the turbu- 
lence quantities. 

Conclusions 

The k-e model of turbulence is used to 
obtain numerical solutions of an axisymmetric turbu- 
lent jet discharging axially into a cylindrical 
tank and directed towards the liquid vapor inter- 
face. The numerical results presented in the paper 
were performed using input conditions representa- 
tive to those of Refs. 10 and 11. In the region 
where the turbulence is approximately isotropic, 
the calculated values of the turbulent velocity com- 
pare well with the experiment. 

The calculated values of the length scale. 

L = k 3/2/ c d 0 no t agree with the conclusion and 
analytical solution of Ref. 10. If the interface 
is far away from the jet, approximately x s /D > 4.5, 
it does not affect the solution up to x/D ^ 3.6. 
However, the boundary conditions for k and e 
imposed at the interface affect the solution in the 
entire region of diffusive turbulence. The results 
obtained by using the boundary conditions (14) for 
large x s /D apprear to be unrealistic. 

References 

1. Labus, T.L., Aydelott, J.C., and Amling, G.E., 
"Zero-Gravity Venting of Three Refrigerants," 
NASA TN D-7480, 1974. 

2. Bullard, B.R., "Liquid Propellant Thermal 
Conditioning System Test Program," 

LMSC-DI59262 , Lockheed Missiles and Space Co., 
Sunnyvale, CA, July 1972, NASA CR-72971. 


3. Robinson, J.W., Giellis, R.T, and Fester, D.A., 
"Long Term Cryogenic Storage Study," 
AFRPL-TR-83-082 , 1984. (Avail. NTIS, 
AD-B089763) . 

4. Aydelott, J.C., Carney, M.J., and Hochstein, 
J.I, "NASA Lewis Research Center Low-Gravity 
Fluid Management Technology Program," AIAA/GNOS 
Paper 85-002, Nov. 1985. (NASA TM-87145). 

5. Merino, F., Blatt, M.K., and Thies, N.C., 
"Filling of Orbital Fluid Management Systems," 
CASD-NAS-78-01 0 , General Dynami cs/Convai r , 

San Diego, CA, July 1978, NASA CR-159404. 

6. Chato, D.J., "Thermodynamic Modeling of the 
No-Vent Fill Methodology for Transferring 
Cryogens in Low Gravity," AIAA-Paper 88-3403, 
July 1988. (NASA TM-100932). 

7. Aydelott, J.C., "Axial Jet Mixing of Ethanol in 
Spherical Containers During Weightlessness," 
NASA TM X-3380 , 1976. 

8. Aydelott, J.C., "Axial Jet Mixing of Ethanol in 
Cylindrical Containers During Weightlessness," 
NASA TP-1487, 1979. 

9. Aydelott, J.C., "Modeling of Space Vehicle 
Propellant Mixing," NASA TP-2107, 1983. 

10. Sonin, A. A., Shimko, M.A., and Chun, J.H., 
"Vapor Condensation Onto a Turbulent Liquid. 1. 
The Steady Condensation Rate as a Function of 
Liquid-Side Turbulence," International Journal 
of Heat and Mass Transfer , Vol . 29, No. 9, 

1986, pp. 1319-1332. 

11. Helmick, M.R., Khoo, B.C., Brown, J.S., and 
Sonin, A. A., "Vapor Condensation Rate at a 
Turbulent Liquid Interface, for Application to 
Cryogenic Hydrogen," AIAA Paper 88-0559, 

Jan. 1988. 

12. Jones, W.P., and Launder, B.E., "The Prediction 
of Laminarization with a Two-Equation Model of 
Turbulence," International Journal of Heat and 
Mass Transfer , Vol. 15, No. 2, Feb. 1972, 

pp. 301-314. 

13. Rod i , W . , Turbulence Models and Their Applica- 
tion i n Hydraul i cs , International Association 
for Hydraulic Research, The Netherlands, 1980, 
pp. 1-92. ' 

14. Launder, B.E., and Spalding, D.B., "The Numeri- 
cal Computation of Turbulent Flows," Comput. 
Methods Appl . Mech. Eng., Vol. 3, No. 2, 

Mar. 1974, pp. 269-289. 

15. Chiappetta, L.M., "User's Manual for a TEACH 
Computer Program for the Analysis of Turbulent, 
Swirling Reacting Flow in a Research Combustor, 
R83-915540-27, United Technologies Research 
Center, East Hartford, CT, Sept. 1983, NASA 
CR-179547. 

16. Escudier, M.P., "The Distribution of Mixing- 
Length in Turbulent Flows Near Walls," Report 
TWF/TN/1 , Imperial College of Science and 
Technology, Dept, of Mechanical Engineering, 
London, 1966. 


5 


ORIGINAL PAGE IS 

02 POOR QUALITY 


17. Syed, S . A . , Chiappetta, L.l 
"Error Reduction Program,'* 
and Whitney Aircraft, East 
Jan. 1985, NASA CR-174775. 


I. , and Gosman, A.D. , 
PWA-5928-25 , Pratt 
Hartford, CT, 


18. Issa. R.I.. "Solution of Implicitly Discretised 
Fluid Flow Equations by Operator-Splitting," 
Journal of Computational Physics , Vol . 62, 

No. 1 , 1986, pp. 40-65. 



FIGURE 1. - GEOMETRIC CONFIGURATION 
AND COORDINATES SYSTEM. 
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FIGURE 2. - COMPARISON OF CALCULATED V VERSUS Q/Dd 
WITH THE EXPERIMENT (REF. 10) AT A FIXED AXIAL 
LOCATION. 



(X s - X)/D 

FIGURE 3. - COMPARISON OF CALCULATED v VERSUS DEPTH 
(X $ - X) BELOW SURFACE WITH EXPERIMENT (REF. 11) AT 
A FIXED JET FLOW RATE. 
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k-e MODEL LENGTH SCALE, L/D in (v) 



AXIAL DISTANCE, x/D 

FIGURE - VARIATION OF v WITH THE AXIAL DISTANCE 
X/D. COMPARISON WITH EXPERIMENT (REF. 10). 



FIGURE 5. - RADIAL DISTRIBUTION OF THE LENGTH SCALE 
L = k 2/3 /E, Rej = 50 000. 
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FIGURE 6. - VARIATION OF TURBULENT KINETIC ENERGY 
WITH THE AXIAL DISTANCE x/D. 



AXIAL DISTANCE. x/D 


FIGURE 7. - VARIATION OF THE LENGTH SCALE, L, WITH THE 
AXIAL DISTANCE x/D. 
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